COVID19 GER

Prepare data

Prevalence

df_ger_prev <- read_csv('Germany_timeseries_prep.csv')

Parsed with column specification:
cols(
  date = col_character(),
  anzahlfall = col_double(),
  kreis = col_double(),
  ewz = col_double(),
  shape__area = col_double(),
  cumcase = col_double(),
  kreis_name = col_character(),
  extra = col_double(),
  agree = col_double(),
  sci = col_double(),
  neuro = col_double(),
  open = col_double(),
  rate_day = col_double(),
  popdens = col_double(),
  runday = col_double()
)

df_ger_prev <- df_ger_prev %>% mutate(date = as.Date(date, "%d%b%Y"),
                                  kreis = as.character(kreis)) %>% 
  dplyr::select(kreis, date, rate_day)

df_ger_prev

Scoial distancing

df_ger_socdist$date %>% summary()

        Min.      1st Qu.       Median         Mean      3rd Qu.         Max. 
"2020-02-25" "2020-03-11" "2020-03-27" "2020-03-27" "2020-04-12" "2020-04-27" 

Personality

df_ger_pers <- read_csv('Germany_timeseries_prep.csv')

Parsed with column specification:
cols(
  date = col_character(),
  anzahlfall = col_double(),
  kreis = col_double(),
  ewz = col_double(),
  shape__area = col_double(),
  cumcase = col_double(),
  kreis_name = col_character(),
  extra = col_double(),
  agree = col_double(),
  sci = col_double(),
  neuro = col_double(),
  open = col_double(),
  rate_day = col_double(),
  popdens = col_double(),
  runday = col_double()
)

df_ger_pers <- df_ger_pers %>% 
  select(kreis, open, sci, extra, agree, neuro) %>%
  dplyr::rename(pers_o = open, 
         pers_c = sci,
         pers_e = extra,
         pers_a = agree,
         pers_n = neuro) %>%
  distinct() %>%
  mutate(kreis = as.character(kreis))

df_ger_pers

NA

Controls

df_ger_ctrl <- read.csv2('Germany_controls.csv', sep = ';', dec=',')

df_ger_ctrl <- df_ger_ctrl %>% select(-kreis_nme) %>%
    mutate(kreis = as.character(kreis),
           popdens = popdens %>% 
             as.character() %>%
             str_replace('\\.', '')%>%
             as.numeric())

df_ger_ctrl

NA

Merge prevalence data

# create sequence of dates
date_sequence <- seq.Date(min(df_ger_prev$date),
                          max(df_ger_prev$date), 1)
                     
# create data frame with time sequence
df_dates = tibble(date_sequence, 1:length(date_sequence)) 
names(df_dates) <- c('date', 'time')


# merge prevalence data
df_ger_prev <- df_ger_prev %>% 
  inner_join(df_ger_pers, by = 'kreis') %>%
  inner_join(df_ger_ctrl, by = 'kreis') %>%
  merge(df_dates, by='date') %>% 
  arrange(kreis)

df_ger_prev

Merge socdist data

# create sequence of dates
date_sequence <- seq.Date(min(df_ger_socdist$date),
                          max(df_ger_socdist$date), 1)
                     
# create data frame with time sequence
df_dates = tibble(date_sequence, 1:length(date_sequence)) 
names(df_dates) <- c('date', 'time')

# merge socdist data
df_ger_socdist <- df_ger_socdist %>% 
  inner_join(df_ger_pers, by = 'kreis') %>%
  inner_join(df_ger_ctrl, by = 'kreis') %>% 
  merge(df_dates, by='date') %>% 
  arrange(kreis)

df_ger_socdist

NA

Control for weekend effect

easter <- seq.Date(as.Date('2020-04-10'), as.Date('2020-04-13'), 1)

df_ger_loess <- df_ger_socdist %>% 
  mutate(weekday = format(date, '%u')) %>% 
  filter(!(weekday %in% c('6','7') | date %in% easter)) %>% 
  split(.$kreis) %>%
  map(~ loess(socdist_single_tile ~ time, data = .)) %>%
  map(predict, 1:max(df_ger_socdist$time)) %>% 
  bind_rows() %>% 
  gather(key = 'kreis', value = 'loess') %>% 
  group_by(kreis) %>% 
  mutate(time = row_number())

df_ger_socdist <- df_ger_socdist %>% merge(df_ger_loess, by=c('kreis', 'time')) %>% 
  mutate(weekday = format(date, '%u')) %>% 
  mutate(socdist_single_tile_clean = ifelse(weekday %in% c('6','7') | date %in% easter, 
                                            loess, socdist_single_tile)) %>%
  arrange(kreis, time) %>% 
  select(-weekday)

df_ger_socdist %>% drop_na()

NA

Explore data

Plot prevalence over time

df_ger_prev %>% ggplot(aes(x=time, y=rate_day)) + 
  geom_point(aes(col=kreis, size=popdens)) + 
  geom_smooth(method="loess", se=T) + 
  theme(legend.position="none") +
  ggtitle("Overall prevalence over time")

pers <- c('pers_o', 'pers_c', 'pers_e', 'pers_a', 'pers_n')

for (i in pers){

gg <- df_ger_prev %>% mutate(prev_tail = cut(.[[i]], 
      breaks = c(-Inf, quantile(.[[i]], 0.2), quantile(.[[i]], 0.8), Inf),
      labels = c('lower tail', 'center', 'upper tail'))) %>% 
  filter(prev_tail != 'center') %>%
  ggplot(aes(x=time, y=rate_day)) + 
  geom_point(aes(col=kreis, size=popdens)) + 
  geom_smooth(method="loess", se=T) + 
  facet_wrap(~prev_tail) + 
  theme(legend.position="none") +
  ggtitle(i)

print(gg)
}

Plot social distancing (single tile) over time

df_ger_socdist %>% ggplot(aes(x=time, y=socdist_single_tile_clean)) + 
  geom_point(aes(col=kreis, size=popdens)) + 
  geom_smooth(method="loess", se=T) + 
  theme(legend.position="none") +
  ggtitle("Overall social distancing (single tile) over time")

pers <- c('pers_o', 'pers_c', 'pers_e', 'pers_a', 'pers_n')

for (i in pers){

gg <- df_ger_socdist %>% mutate(prev_tail = cut(.[[i]], 
      breaks = c(-Inf, quantile(.[[i]], 0.2), quantile(.[[i]], 0.8), Inf),
      labels = c('lower tail', 'center', 'upper tail'))) %>% 
  filter(prev_tail != 'center') %>%
  ggplot(aes(x=time, y=socdist_single_tile_clean)) + 
  geom_point(aes(col=kreis, size=popdens)) + 
  geom_smooth(method="loess", se=T) + 
  facet_wrap(~prev_tail) + 
  theme(legend.position="none") +
  ggtitle(i)

print(gg)
}

df_ger_socdist <- df_ger_socdist %>% mutate(socdist_single_tile = socdist_single_tile_clean) %>% 
  select(-loess, -socdist_single_tile_clean)

Correlations

df_ger_prev %>% group_by(kreis) %>% 
  summarize_if(is.numeric, mean) %>% 
  select(-kreis, -time) %>% 
  cor(use='pairwise.complete') %>% 
  round(3) %>% as.data.frame()

df_ger_socdist %>% group_by(kreis) %>% 
  summarize_if(is.numeric, mean) %>% 
  select(-kreis, -time) %>% 
  cor(use='pairwise.complete') %>% 
  round(3) %>% as.data.frame()

NA
NA

Rescale data

lvl2_scaled <- df_ger_prev %>% 
  dplyr::select(-time, -rate_day, -date) %>% 
  distinct() %>% 
  mutate_at(vars(-kreis), scale)
  
lvl1_scaled <- df_ger_prev %>% select(kreis, time, rate_day)

df_ger_prev_scaled <- plyr::join(lvl1_scaled, lvl2_scaled, by = 'kreis')

df_ger_prev_scaled

NA

lvl2_scaled <- df_ger_socdist %>% 
  dplyr::select(-time, -socdist_single_tile, -date) %>% 
  distinct() %>% 
  mutate_at(vars(-kreis), scale)
  
lvl1_scaled <- df_ger_socdist %>% select(kreis, time, socdist_single_tile)

df_ger_socdist_scaled <- plyr::join(lvl1_scaled, lvl2_scaled, by = 'kreis')

df_ger_socdist_scaled

NA

Predict Prevalence

Extract first day of covid outbreak

# get onset day
df_ger_onset_prev <- df_ger_prev_scaled %>% 
  group_by(kreis) %>% 
  mutate(rate_cs = cumsum(rate_day)) %>% 
  filter(rate_cs > 0) %>%
  summarize(onset_prev = min(time))
  
# merge with county data
df_ger_onset_prev <- df_ger_prev_scaled %>% 
  select(-time, -rate_day) %>%
  distinct() %>% 
  left_join(df_ger_onset_prev, by = 'kreis')

# handle censored data
df_ger_onset_prev <- df_ger_onset_prev %>% 
  mutate(event = ifelse(is.na(onset_prev), 0, 1)) %>% 
  mutate(onset_prev = replace_na(onset_prev, as.numeric(diff(range(df_ger_prev$date)))+1))

df_ger_onset_prev

NA

Extract slopes

# cut time series before onset
df_ger_prev_scaled <- df_ger_prev_scaled %>% 
  group_by(kreis) %>% 
  mutate(rate_cs = cumsum(rate_day)) %>% 
  filter(rate_cs > 0) %>%
  mutate(time = time-min(time)+1) %>%
  ungroup() %>%
  filter(time <= 30) %>%
  select(-rate_cs)

# drop counties with little data
df_ger_prev_scaled <- df_ger_prev_scaled %>%
  group_by(kreis) %>%
  filter(n() == 30) %>%
  ungroup()

# extract slope prevalence
df_ger_slope_prev <- df_ger_prev_scaled %>% 
  split(.$kreis) %>% 
  map(~ lm(rate_day ~ time, data = .)) %>%
  map(coef) %>% 
  map_dbl('time') %>% 
  as.data.frame() %>% 
  rownames_to_column('kreis') %>% 
  rename(slope_prev = '.')

# merge with county data
df_ger_slope_prev <- df_ger_prev_scaled %>% 
  select(-time, -rate_day) %>%
  distinct() %>% 
  inner_join(df_ger_slope_prev, by = 'kreis') %>%
  drop_na()

# standardize slopes
df_ger_slope_prev <- df_ger_slope_prev %>% 
  mutate(slope_prev = scale(slope_prev))

Explore distributions

df_ger_onset_prev %>% ggplot(aes(onset_prev)) + geom_histogram()

df_ger_slope_prev %>% ggplot(aes(slope_prev)) + geom_histogram()

Predict COVID onset with time-to-event regression

# predict onset from personality
cox_onset_prev <- coxph(Surv(onset_prev, event) ~ 
                          pers_o + pers_c + pers_e + pers_a + pers_n, 
                        data = df_ger_onset_prev)
cox_onset_prev %>% summary()

Call:
coxph(formula = Surv(onset_prev, event) ~ pers_o + pers_c + pers_e + 
    pers_a + pers_n, data = df_ger_onset_prev)

  n= 400, number of events= 400 

            coef exp(coef)  se(coef)      z Pr(>|z|)    
pers_o  0.131653  1.140712  0.051953  2.534   0.0113 *  
pers_c -0.104500  0.900775  0.053314 -1.960   0.0500 *  
pers_e  0.056186  1.057794  0.052903  1.062   0.2882    
pers_a -0.008624  0.991413  0.049615 -0.174   0.8620    
pers_n -0.276799  0.758207  0.057728 -4.795 1.63e-06 ***
---
Signif. codes:  0 ‘***’ 0.001 ‘**’ 0.01 ‘*’ 0.05 ‘.’ 0.1 ‘ ’ 1

       exp(coef) exp(-coef) lower .95 upper .95
pers_o    1.1407     0.8766    1.0303     1.263
pers_c    0.9008     1.1102    0.8114     1.000
pers_e    1.0578     0.9454    0.9536     1.173
pers_a    0.9914     1.0087    0.8995     1.093
pers_n    0.7582     1.3189    0.6771     0.849

Concordance= 0.592  (se = 0.018 )
Likelihood ratio test= 39.92  on 5 df,   p=2e-07
Wald test            = 37.99  on 5 df,   p=4e-07
Score (logrank) test = 37.99  on 5 df,   p=4e-07

# predict onset from personality with controls
cox_onset_prev_ctrl <- coxph(Surv(onset_prev, event) ~ 
                               pers_o + pers_c + pers_e + pers_a + pers_n +
                               women + academics + afd + hospital_beds +
                               tourism_beds + gdp + manufact + airport + age +
                               popdens,
                             data = df_ger_onset_prev)
cox_onset_prev_ctrl %>% summary()

Call:
coxph(formula = Surv(onset_prev, event) ~ pers_o + pers_c + pers_e + 
    pers_a + pers_n + women + academics + afd + hospital_beds + 
    tourism_beds + gdp + manufact + airport + age + popdens, 
    data = df_ger_onset_prev)

  n= 392, number of events= 392 
   (8 observations deleted due to missingness)

                   coef exp(coef)  se(coef)      z Pr(>|z|)    
pers_o         0.050510  1.051807  0.059853  0.844 0.398726    
pers_c        -0.087556  0.916168  0.061192 -1.431 0.152473    
pers_e         0.038883  1.039649  0.057895  0.672 0.501825    
pers_a         0.008022  1.008054  0.057017  0.141 0.888110    
pers_n        -0.222357  0.800630  0.063469 -3.503 0.000459 ***
women          0.066190  1.068429  0.067498  0.981 0.326785    
academics      0.267424  1.306595  0.094669  2.825 0.004730 ** 
afd           -0.072815  0.929773  0.075881 -0.960 0.337260    
hospital_beds -0.236175  0.789642  0.069836 -3.382 0.000720 ***
tourism_beds   0.076981  1.080022  0.054241  1.419 0.155829    
gdp           -0.175570  0.838979  0.115831 -1.516 0.129586    
manufact       0.147898  1.159395  0.094458  1.566 0.117404    
airport       -0.171432  0.842457  0.061687 -2.779 0.005451 ** 
age           -0.176055  0.838572  0.085439 -2.061 0.039343 *  
popdens        0.062986  1.065012  0.082801  0.761 0.446841    
---
Signif. codes:  0 ‘***’ 0.001 ‘**’ 0.01 ‘*’ 0.05 ‘.’ 0.1 ‘ ’ 1

              exp(coef) exp(-coef) lower .95 upper .95
pers_o           1.0518     0.9507    0.9354    1.1827
pers_c           0.9162     1.0915    0.8126    1.0329
pers_e           1.0396     0.9619    0.9281    1.1646
pers_a           1.0081     0.9920    0.9015    1.1272
pers_n           0.8006     1.2490    0.7070    0.9067
women            1.0684     0.9360    0.9360    1.2196
academics        1.3066     0.7653    1.0853    1.5730
afd              0.9298     1.0755    0.8013    1.0789
hospital_beds    0.7896     1.2664    0.6886    0.9055
tourism_beds     1.0800     0.9259    0.9711    1.2012
gdp              0.8390     1.1919    0.6686    1.0528
manufact         1.1594     0.8625    0.9634    1.3952
airport          0.8425     1.1870    0.7465    0.9507
age              0.8386     1.1925    0.7093    0.9914
popdens          1.0650     0.9390    0.9055    1.2527

Concordance= 0.658  (se = 0.016 )
Likelihood ratio test= 101.3  on 15 df,   p=7e-15
Wald test            = 98.22  on 15 df,   p=3e-14
Score (logrank) test = 102.1  on 15 df,   p=5e-15

Predict prevalence slopes with linear models

# predict slopes from personality
lm_slope_prev <- lm(slope_prev ~ pers_o + pers_c + pers_e + pers_a + pers_n, 
                         data = df_ger_slope_prev)
lm_slope_prev %>% summary()

Call:
lm(formula = slope_prev ~ pers_o + pers_c + pers_e + pers_a + 
    pers_n, data = df_ger_slope_prev)

Residuals:
    Min      1Q  Median      3Q     Max 
-1.1994 -0.5725 -0.1888  0.2940  8.6104 

Coefficients:
            Estimate Std. Error t value Pr(>|t|)
(Intercept)  0.00174    0.05054   0.034    0.973
pers_o      -0.03181    0.05420  -0.587    0.558
pers_c      -0.05989    0.05799  -1.033    0.302
pers_e       0.07814    0.05745   1.360    0.175
pers_a       0.09131    0.05839   1.564    0.119
pers_n       0.06209    0.05971   1.040    0.299

Residual standard error: 1 on 386 degrees of freedom
Multiple R-squared:  0.01215,   Adjusted R-squared:  -0.0006472 
F-statistic: 0.9494 on 5 and 386 DF,  p-value: 0.4489

lm_slope_prev %>% confint(level=0.9)

                     5 %       95 %
(Intercept) -0.081594195 0.08507360
pers_o      -0.121173276 0.05756145
pers_c      -0.155494307 0.03572519
pers_e      -0.016586494 0.17286042
pers_a      -0.004968057 0.18758968
pers_n      -0.036363942 0.16053810

# predict slopes from personality with controls
lm_slope_prev_ctrl <- lm(slope_prev ~  
                               pers_o + pers_c + pers_e + pers_a + pers_n +
                               women + academics + afd + hospital_beds +
                               tourism_beds + gdp + manufact + airport + age +
                               popdens,
                         data = df_ger_slope_prev)
lm_slope_prev_ctrl %>% summary()

Call:
lm(formula = slope_prev ~ pers_o + pers_c + pers_e + pers_a + 
    pers_n + women + academics + afd + hospital_beds + tourism_beds + 
    gdp + manufact + airport + age + popdens, data = df_ger_slope_prev)

Residuals:
    Min      1Q  Median      3Q     Max 
-1.4465 -0.5788 -0.0848  0.3212  8.0256 

Coefficients:
                Estimate Std. Error t value Pr(>|t|)  
(Intercept)   -0.0005982  0.0488695  -0.012   0.9902  
pers_o        -0.0160080  0.0603371  -0.265   0.7909  
pers_c         0.0018010  0.0584550   0.031   0.9754  
pers_e         0.0769996  0.0572311   1.345   0.1793  
pers_a         0.0832263  0.0586703   1.419   0.1569  
pers_n         0.0621839  0.0593346   1.048   0.2953  
women          0.1316729  0.0715869   1.839   0.0667 .
academics     -0.0618388  0.0955468  -0.647   0.5179  
afd           -0.0225847  0.0712477  -0.317   0.7514  
hospital_beds  0.0939475  0.0657371   1.429   0.1538  
tourism_beds  -0.0400654  0.0550991  -0.727   0.4676  
gdp           -0.0130746  0.1126142  -0.116   0.9076  
manufact       0.1729185  0.0894459   1.933   0.0540 .
airport        0.1503036  0.0612082   2.456   0.0145 *
age           -0.1912130  0.0908032  -2.106   0.0359 *
popdens       -0.1393613  0.0790525  -1.763   0.0787 .
---
Signif. codes:  0 ‘***’ 0.001 ‘**’ 0.01 ‘*’ 0.05 ‘.’ 0.1 ‘ ’ 1

Residual standard error: 0.9665 on 376 degrees of freedom
Multiple R-squared:  0.1016,    Adjusted R-squared:  0.06581 
F-statistic: 2.836 on 15 and 376 DF,  p-value: 0.0003208

lm_slope_prev_ctrl %>% confint(level=0.9)

                      5 %         95 %
(Intercept)   -0.08118002  0.079983552
pers_o        -0.11549884  0.083482831
pers_c        -0.09458638  0.098188366
pers_e        -0.01736960  0.171368869
pers_a        -0.01351600  0.179968683
pers_n        -0.03565397  0.160021700
women          0.01363199  0.249713732
academics     -0.21938756  0.095709935
afd           -0.14006616  0.094896812
hospital_beds -0.01444736  0.202342442
tourism_beds  -0.13091924  0.050788398
gdp           -0.19876601  0.172616820
manufact       0.02542968  0.320407240
airport        0.04937648  0.251230766
age           -0.34093993 -0.041486083
popdens       -0.26971227 -0.009010392

CRF predicting slopes

ctrls <- cforest_unbiased(ntree=500, mtry=5)

crf_slope_prev <- cforest(slope_prev ~  
                               pers_o + pers_c + pers_e + pers_a + pers_n +
                               women + academics + afd + hospital_beds +
                               tourism_beds + gdp + manufact + airport + age +
                               popdens,
                           data = df_ger_slope_prev, 
                         controls = ctrls)

crf_slope_prev_varimp <- varimp(crf_slope_prev, nperm = 1)
crf_slope_prev_varimp_cond <- varimp(crf_slope_prev, conditional = T, nperm = 1)

crf_slope_prev_varimp

       pers_o        pers_c        pers_e        pers_a        pers_n         women     academics 
  0.001488878  -0.004182147   0.001303730   0.012206527   0.001014518   0.009796793   0.003829958 
          afd hospital_beds  tourism_beds           gdp      manufact       airport           age 
  0.015500844   0.018005690  -0.001781128   0.021456347   0.035828919   0.017943409   0.023111995 
      popdens 
  0.011791481 

crf_slope_prev_varimp %>% as.data.frame() %>% 
  rownames_to_column('variable') %>%
  ggplot(aes(x=variable, y=.)) +
  geom_bar(stat = 'identity') + 
  theme(axis.text.x = element_text(angle = 90))

crf_slope_prev_varimp_cond

       pers_o        pers_c        pers_e        pers_a        pers_n         women     academics 
-0.0002143011 -0.0014382430 -0.0047726834  0.0097148040  0.0009568305  0.0163477254  0.0051396731 
          afd hospital_beds  tourism_beds           gdp      manufact       airport           age 
 0.0163462437  0.0164595151 -0.0030080531  0.0227972666  0.0308178428  0.0157462891  0.0238096354 
      popdens 
 0.0150246388 

crf_slope_prev_varimp_cond %>% as.data.frame() %>% 
  rownames_to_column('variable') %>%
  ggplot(aes(x=variable, y=.)) +
  geom_bar(stat = 'identity') +
  theme(axis.text.x = element_text(angle = 90))

Predict Social Distancing

Change point analysis

# keep only counties with full data
kreis_complete <- df_ger_socdist_scaled %>% 
  group_by(kreis) %>% 
  summarize(n = n()) %>% 
  filter(n==max(.$n)) %>% 
  .$kreis

# run changepoint analysis
df_ger_socdist_cpt_results <- df_ger_socdist_scaled %>% 
  select(kreis, socdist_single_tile) %>%
  filter(kreis %in% kreis_complete) %>% 
  split(.$kreis) %>%
  map(~ cpt.meanvar(as.vector(.$socdist_single_tile),
                    #penalty = 'Asymptotic',
                    class=TRUE,
                    param.estimates=TRUE,
                    Q=1,
                    test.stat = 'Normal'))

# calculate change point
df_ger_socdist_cpt_day <- df_ger_socdist_cpt_results %>% 
  map(cpts) %>% 
  unlist() %>% 
  as.data.frame() %>% 
  rename(cpt_day_socdist = '.') %>%
  rownames_to_column('kreis')

# calculate mean differences
df_ger_socdist_cpt_mean_diff <- df_ger_socdist_cpt_results %>% 
  map(param.est) %>% 
  map(~ .$mean) %>% 
  map(~ .[2]-.[1]) %>% 
  unlist() %>% 
  as.data.frame() %>% 
  rename(mean_diff_socdist = '.') %>%
  rownames_to_column('kreis')

# calculate varaince differences
df_ger_socdist_cpt_var_diff <- df_ger_socdist_cpt_results %>% 
  map(param.est) %>% 
  map(~ .$variance) %>% 
  map(~ .[2]-.[1]) %>% 
  unlist() %>% 
  as.data.frame() %>% 
  rename(var_diff_socdist = '.') %>%
  rownames_to_column('kreis')

# merge with county data
df_ger_cpt_socdist <- df_ger_socdist_scaled %>% 
  select(-time, -socdist_single_tile) %>%
  distinct() %>% 
  left_join(df_ger_socdist_cpt_day, by='kreis') %>%
  left_join(df_ger_socdist_cpt_mean_diff, by='kreis') %>%
  left_join(df_ger_socdist_cpt_var_diff, by='kreis') %>%
  left_join(select(df_ger_onset_prev, kreis, onset_prev), by='kreis') %>%
  left_join(select(df_ger_slope_prev, kreis, slope_prev), by='kreis') 

# standardize mean/var differences
df_ger_cpt_socdist <- df_ger_cpt_socdist %>% 
  mutate(mean_diff_socdist = scale(mean_diff_socdist),
         var_diff_socdist = scale(var_diff_socdist))

# handle censored data
df_ger_cpt_socdist <- df_ger_cpt_socdist %>% 
  mutate(cpt_day_socdist = ifelse(is.na(cpt_day_socdist), 
                                  as.numeric(diff(range(df_ger_socdist$date))), 
                                  cpt_day_socdist)) %>% 
  mutate(event = ifelse(cpt_day_socdist >= 60, 0, 1))

df_ger_cpt_socdist$cpt_day_socdist %>% hist()

df_ger_cpt_socdist$mean_diff_socdist %>% hist()

df_ger_cpt_socdist$var_diff_socdist %>% hist()

for(i in head(df_ger_socdist_cpt_results, 5)){
  plot(i)
}

NA

Predicting change points with time-to-event regression

# predict hazard from personality
cox_cpt_socdist <- coxph(Surv(cpt_day_socdist, event) ~ 
                           pers_o + pers_c + pers_e + pers_a + pers_n, 
                  data = df_ger_cpt_socdist)
cox_cpt_socdist %>% summary()

Call:
coxph(formula = Surv(cpt_day_socdist, event) ~ pers_o + pers_c + 
    pers_e + pers_a + pers_n, data = df_ger_cpt_socdist)

  n= 400, number of events= 400 

           coef exp(coef) se(coef)      z Pr(>|z|)  
pers_o  0.13238   1.14154  0.05555  2.383   0.0172 *
pers_c -0.04634   0.95472  0.05364 -0.864   0.3876  
pers_e -0.08928   0.91459  0.05443 -1.640   0.1009  
pers_a -0.05236   0.94898  0.04881 -1.073   0.2833  
pers_n -0.08884   0.91499  0.05376 -1.653   0.0984 .
---
Signif. codes:  0 ‘***’ 0.001 ‘**’ 0.01 ‘*’ 0.05 ‘.’ 0.1 ‘ ’ 1

       exp(coef) exp(-coef) lower .95 upper .95
pers_o    1.1415      0.876    1.0238     1.273
pers_c    0.9547      1.047    0.8594     1.061
pers_e    0.9146      1.093    0.8221     1.018
pers_a    0.9490      1.054    0.8624     1.044
pers_n    0.9150      1.093    0.8235     1.017

Concordance= 0.595  (se = 0.024 )
Likelihood ratio test= 10.81  on 5 df,   p=0.06
Wald test            = 10.76  on 5 df,   p=0.06
Score (logrank) test = 10.71  on 5 df,   p=0.06

# predict hazard from personality with controls
cox_cpt_socdist_ctrl <- coxph(Surv(cpt_day_socdist, event) ~ 
                                pers_o + pers_c + pers_e + pers_a + pers_n +
                               women + academics + afd + hospital_beds +
                               tourism_beds + gdp + manufact + airport + age +
                               popdens,
                           data = df_ger_cpt_socdist,)
cox_cpt_socdist_ctrl %>% summary()

Call:
coxph(formula = Surv(cpt_day_socdist, event) ~ pers_o + pers_c + 
    pers_e + pers_a + pers_n + women + academics + afd + hospital_beds + 
    tourism_beds + gdp + manufact + airport + age + popdens, 
    data = df_ger_cpt_socdist)

  n= 392, number of events= 392 
   (8 observations deleted due to missingness)

                  coef exp(coef) se(coef)      z Pr(>|z|)   
pers_o         0.01147   1.01154  0.06432  0.178  0.85846   
pers_c        -0.04040   0.96041  0.06122 -0.660  0.50935   
pers_e        -0.16691   0.84628  0.06130 -2.723  0.00647 **
pers_a        -0.02162   0.97861  0.05627 -0.384  0.70084   
pers_n        -0.05115   0.95013  0.06007 -0.852  0.39442   
women          0.10582   1.11163  0.07640  1.385  0.16603   
academics     -0.01816   0.98201  0.10793 -0.168  0.86640   
afd           -0.03782   0.96288  0.08037 -0.471  0.63792   
hospital_beds -0.19535   0.82255  0.06587 -2.965  0.00302 **
tourism_beds   0.08410   1.08774  0.05639  1.492  0.13583   
gdp            0.14418   1.15509  0.12603  1.144  0.25262   
manufact      -0.14887   0.86168  0.09244 -1.610  0.10731   
airport       -0.09012   0.91382  0.06467 -1.394  0.16343   
age           -0.14891   0.86165  0.09656 -1.542  0.12304   
popdens        0.17167   1.18729  0.08440  2.034  0.04196 * 
---
Signif. codes:  0 ‘***’ 0.001 ‘**’ 0.01 ‘*’ 0.05 ‘.’ 0.1 ‘ ’ 1

              exp(coef) exp(-coef) lower .95 upper .95
pers_o           1.0115     0.9886    0.8917    1.1474
pers_c           0.9604     1.0412    0.8518    1.0829
pers_e           0.8463     1.1816    0.7505    0.9543
pers_a           0.9786     1.0219    0.8764    1.0927
pers_n           0.9501     1.0525    0.8446    1.0688
women            1.1116     0.8996    0.9570    1.2912
academics        0.9820     1.0183    0.7948    1.2134
afd              0.9629     1.0385    0.8225    1.1272
hospital_beds    0.8226     1.2157    0.7229    0.9359
tourism_beds     1.0877     0.9193    0.9739    1.2148
gdp              1.1551     0.8657    0.9023    1.4787
manufact         0.8617     1.1605    0.7189    1.0328
airport          0.9138     1.0943    0.8050    1.0373
age              0.8616     1.1606    0.7131    1.0412
popdens          1.1873     0.8423    1.0063    1.4009

Concordance= 0.697  (se = 0.022 )
Likelihood ratio test= 62.2  on 15 df,   p=1e-07
Wald test            = 66.91  on 15 df,   p=2e-08
Score (logrank) test = 68.7  on 15 df,   p=8e-09

# predict hazard from personality with controls
cox_cpt_socdist_ctrl2 <- coxph(Surv(cpt_day_socdist, event) ~ 
                                 pers_o + pers_c + pers_e + pers_a + pers_n +
                               women + academics + afd + hospital_beds +
                               tourism_beds + gdp + manufact + airport + age +
                               popdens + onset_prev + slope_prev,
                  data = df_ger_cpt_socdist)
cox_cpt_socdist_ctrl2 %>% summary()

Call:
coxph(formula = Surv(cpt_day_socdist, event) ~ pers_o + pers_c + 
    pers_e + pers_a + pers_n + women + academics + afd + hospital_beds + 
    tourism_beds + gdp + manufact + airport + age + popdens + 
    onset_prev + slope_prev, data = df_ger_cpt_socdist)

  n= 392, number of events= 392 
   (8 observations deleted due to missingness)

                   coef exp(coef)  se(coef)      z Pr(>|z|)   
pers_o         0.009217  1.009259  0.064871  0.142  0.88702   
pers_c        -0.041231  0.959607  0.061433 -0.671  0.50212   
pers_e        -0.172081  0.841911  0.062033 -2.774  0.00554 **
pers_a        -0.014363  0.985739  0.057276 -0.251  0.80199   
pers_n        -0.037497  0.963197  0.061112 -0.614  0.53949   
women          0.121990  1.129743  0.077951  1.565  0.11759   
academics     -0.041589  0.959264  0.109525 -0.380  0.70416   
afd           -0.021631  0.978601  0.081267 -0.266  0.79010   
hospital_beds -0.189670  0.827232  0.066132 -2.868  0.00413 **
tourism_beds   0.082883  1.086415  0.056235  1.474  0.14052   
gdp            0.175596  1.191956  0.127605  1.376  0.16879   
manufact      -0.170477  0.843263  0.093993 -1.814  0.06972 . 
airport       -0.077842  0.925111  0.065550 -1.188  0.23502   
age           -0.149096  0.861486  0.098774 -1.509  0.13118   
popdens        0.157912  1.171064  0.084451  1.870  0.06150 . 
onset_prev    -0.006740  0.993283  0.005197 -1.297  0.19468   
slope_prev    -0.026992  0.973369  0.064662 -0.417  0.67637   
---
Signif. codes:  0 ‘***’ 0.001 ‘**’ 0.01 ‘*’ 0.05 ‘.’ 0.1 ‘ ’ 1

              exp(coef) exp(-coef) lower .95 upper .95
pers_o           1.0093     0.9908    0.8888    1.1461
pers_c           0.9596     1.0421    0.8507    1.0824
pers_e           0.8419     1.1878    0.7455    0.9508
pers_a           0.9857     1.0145    0.8811    1.1028
pers_n           0.9632     1.0382    0.8545    1.0858
women            1.1297     0.8852    0.9697    1.3162
academics        0.9593     1.0425    0.7739    1.1890
afd              0.9786     1.0219    0.8345    1.1476
hospital_beds    0.8272     1.2089    0.7267    0.9417
tourism_beds     1.0864     0.9205    0.9730    1.2130
gdp              1.1920     0.8390    0.9282    1.5307
manufact         0.8433     1.1859    0.7014    1.0138
airport          0.9251     1.0810    0.8136    1.0519
age              0.8615     1.1608    0.7099    1.0455
popdens          1.1711     0.8539    0.9924    1.3819
onset_prev       0.9933     1.0068    0.9832    1.0035
slope_prev       0.9734     1.0274    0.8575    1.1049

Concordance= 0.7  (se = 0.022 )
Likelihood ratio test= 65.24  on 17 df,   p=1e-07
Wald test            = 69.97  on 17 df,   p=2e-08
Score (logrank) test = 71.77  on 17 df,   p=1e-08

Linear models predicting mean differences

lm_meandiff_socdist <- lm(mean_diff_socdist ~ 
                            pers_o + pers_c + pers_e + pers_a + pers_n, 
                         data = df_ger_cpt_socdist)
lm_meandiff_socdist %>% summary()

Call:
lm(formula = mean_diff_socdist ~ pers_o + pers_c + pers_e + pers_a + 
    pers_n, data = df_ger_cpt_socdist)

Residuals:
    Min      1Q  Median      3Q     Max 
-2.6704 -0.6077 -0.0112  0.5246  3.9074 

Coefficients:
              Estimate Std. Error t value Pr(>|t|)    
(Intercept) -5.881e-16  4.603e-02   0.000   1.0000    
pers_o       3.573e-01  4.928e-02   7.251 2.21e-12 ***
pers_c      -5.383e-02  5.227e-02  -1.030   0.3037    
pers_e       5.436e-02  5.253e-02   1.035   0.3014    
pers_a      -1.307e-02  5.224e-02  -0.250   0.8026    
pers_n      -1.268e-01  5.423e-02  -2.338   0.0199 *  
---
Signif. codes:  0 ‘***’ 0.001 ‘**’ 0.01 ‘*’ 0.05 ‘.’ 0.1 ‘ ’ 1

Residual standard error: 0.9206 on 394 degrees of freedom
Multiple R-squared:  0.1631,    Adjusted R-squared:  0.1525 
F-statistic: 15.36 on 5 and 394 DF,  p-value: 8.43e-14

lm_meandiff_socdist %>% confint(level=0.9)

                    5 %        95 %
(Intercept) -0.07589246  0.07589246
pers_o       0.27605447  0.43854018
pers_c      -0.14000439  0.03234850
pers_e      -0.03224219  0.14095947
pers_a      -0.09919902  0.07306852
pers_n      -0.21617769 -0.03736628

lm_meandiff_socdist_ctrl <- lm(mean_diff_socdist ~ 
                                 pers_o + pers_c + pers_e + pers_a + pers_n +
                               women + academics + afd + hospital_beds +
                               tourism_beds + gdp + manufact + airport + age +
                               popdens,
                           data = df_ger_cpt_socdist,)
lm_meandiff_socdist_ctrl %>% summary()

Call:
lm(formula = mean_diff_socdist ~ pers_o + pers_c + pers_e + pers_a + 
    pers_n + women + academics + afd + hospital_beds + tourism_beds + 
    gdp + manufact + airport + age + popdens, data = df_ger_cpt_socdist)

Residuals:
     Min       1Q   Median       3Q      Max 
-2.82186 -0.48179 -0.06099  0.41463  2.54031 

Coefficients:
               Estimate Std. Error t value Pr(>|t|)    
(Intercept)   -0.005219   0.036964  -0.141 0.887799    
pers_o         0.106164   0.045638   2.326 0.020538 *  
pers_c         0.048613   0.044214   1.099 0.272256    
pers_e        -0.037110   0.043289  -0.857 0.391844    
pers_a        -0.004674   0.044377  -0.105 0.916179    
pers_n         0.014518   0.044880   0.323 0.746510    
women          0.124196   0.054147   2.294 0.022360 *  
academics      0.081958   0.072270   1.134 0.257493    
afd           -0.146522   0.053891  -2.719 0.006854 ** 
hospital_beds -0.179132   0.049722  -3.603 0.000357 ***
tourism_beds   0.026876   0.041676   0.645 0.519393    
gdp            0.322525   0.085179   3.786 0.000178 ***
manufact      -0.044887   0.067655  -0.663 0.507442    
airport       -0.059031   0.046297  -1.275 0.203079    
age           -0.326924   0.068682  -4.760 2.77e-06 ***
popdens       -0.016479   0.059794  -0.276 0.783013    
---
Signif. codes:  0 ‘***’ 0.001 ‘**’ 0.01 ‘*’ 0.05 ‘.’ 0.1 ‘ ’ 1

Residual standard error: 0.7311 on 376 degrees of freedom
  (8 observations deleted due to missingness)
Multiple R-squared:  0.4871,    Adjusted R-squared:  0.4666 
F-statistic: 23.81 on 15 and 376 DF,  p-value: < 2.2e-16

lm_meandiff_socdist_ctrl %>% confint(level=0.9)

                      5 %        95 %
(Intercept)   -0.06616945  0.05573185
pers_o         0.03091038  0.18141662
pers_c        -0.02429231  0.12151912
pers_e        -0.10848908  0.03426939
pers_a        -0.07784795  0.06850047
pers_n        -0.05948507  0.08852057
women          0.03491147  0.21347956
academics     -0.03720884  0.20112538
afd           -0.23538302 -0.05766115
hospital_beds -0.26111965 -0.09714365
tourism_beds  -0.04184386  0.09559660
gdp            0.18207149  0.46297894
manufact      -0.15644458  0.06667126
airport       -0.13537018  0.01730886
age           -0.44017449 -0.21367285
popdens       -0.11507384  0.08211648

lm_meandiff_socdist_ctrl2 <- lm(mean_diff_socdist ~ 
                                  pers_o + pers_c + pers_e + pers_a + pers_n +
                               women + academics + afd + hospital_beds +
                               tourism_beds + gdp + manufact + airport + age +
                               popdens + onset_prev + slope_prev,
                  data = df_ger_cpt_socdist)
lm_meandiff_socdist_ctrl2 %>% summary()

Call:
lm(formula = mean_diff_socdist ~ pers_o + pers_c + pers_e + pers_a + 
    pers_n + women + academics + afd + hospital_beds + tourism_beds + 
    gdp + manufact + airport + age + popdens + onset_prev + slope_prev, 
    data = df_ger_cpt_socdist)

Residuals:
    Min      1Q  Median      3Q     Max 
-2.7924 -0.4819 -0.0360  0.4080  2.6169 

Coefficients:
               Estimate Std. Error t value Pr(>|t|)    
(Intercept)    0.264112   0.201982   1.308 0.191811    
pers_o         0.108998   0.044269   2.462 0.014259 *  
pers_c         0.051255   0.042940   1.194 0.233373    
pers_e        -0.056367   0.042250  -1.334 0.182976    
pers_a        -0.020230   0.043157  -0.469 0.639515    
pers_n         0.009412   0.043847   0.215 0.830160    
women          0.098099   0.052769   1.859 0.063809 .  
academics      0.084990   0.070407   1.207 0.228148    
afd           -0.135862   0.052495  -2.588 0.010026 *  
hospital_beds -0.190103   0.048586  -3.913 0.000108 ***
tourism_beds   0.036003   0.040470   0.890 0.374241    
gdp            0.334790   0.082937   4.037 6.58e-05 ***
manufact      -0.086818   0.066355  -1.308 0.191540    
airport       -0.087019   0.045262  -1.923 0.055297 .  
age           -0.278359   0.067715  -4.111 4.85e-05 ***
popdens        0.009521   0.058234   0.163 0.870216    
onset_prev    -0.004807   0.003549  -1.354 0.176427    
slope_prev     0.210948   0.042472   4.967 1.04e-06 ***
---
Signif. codes:  0 ‘***’ 0.001 ‘**’ 0.01 ‘*’ 0.05 ‘.’ 0.1 ‘ ’ 1

Residual standard error: 0.7091 on 374 degrees of freedom
  (8 observations deleted due to missingness)
Multiple R-squared:  0.5201,    Adjusted R-squared:  0.4983 
F-statistic: 23.84 on 17 and 374 DF,  p-value: < 2.2e-16

lm_meandiff_socdist_ctrl2 %>% confint(level=0.9)

                      5 %         95 %
(Intercept)   -0.06894364  0.597167453
pers_o         0.03600153  0.181994243
pers_c        -0.01955027  0.122060396
pers_e        -0.12603448  0.013301070
pers_a        -0.09139314  0.050932815
pers_n        -0.06288910  0.081712228
women          0.01108590  0.185112674
academics     -0.03110714  0.201088040
afd           -0.22242240 -0.049301807
hospital_beds -0.27021921 -0.109987126
tourism_beds  -0.03072965  0.102735952
gdp            0.19803273  0.471548214
manufact      -0.19623295  0.022596103
airport       -0.16165358 -0.012383791
age           -0.39001650 -0.166702289
popdens       -0.08650307  0.105545112
onset_prev    -0.01065842  0.001045275
slope_prev     0.14091460  0.280982005

CRF predicting mean difference

ctrls <- cforest_unbiased(ntree=500, mtry=5)

crf_meandiff_socdist <- cforest(mean_diff_socdist ~ 
                                  pers_o + pers_c + pers_e + pers_a + pers_n +
                               women + academics + afd + hospital_beds +
                               tourism_beds + gdp + manufact + airport + age +
                               popdens + onset_prev + slope_prev,
                              data = df_ger_cpt_socdist %>% drop_na(),
                         controls = ctrls)

crf_meandiff_socdist_varimp <- varimp(crf_meandiff_socdist, nperm = 1)
crf_meandiff_socdist_varimp_cond <- varimp(crf_meandiff_socdist, conditional = T, nperm = 1)

crf_meandiff_socdist_varimp

       pers_o        pers_c        pers_e        pers_a        pers_n         women     academics 
 0.0180211122  0.0015870098  0.0027043254  0.0005814055  0.0026234653  0.0099984261  0.1482181959 
          afd hospital_beds  tourism_beds           gdp      manufact       airport           age 
 0.1397866414  0.0112907121  0.0004702307  0.1025281549  0.0089789021  0.0165751547  0.2422748719 
      popdens    onset_prev    slope_prev 
 0.0433761825  0.0016622854  0.0148192230 

crf_meandiff_socdist_varimp %>% as.data.frame() %>% 
  rownames_to_column('variable') %>%
  ggplot(aes(x=variable, y=.)) +
  geom_bar(stat = 'identity') +
  theme(axis.text.x = element_text(angle = 90))

crf_meandiff_socdist_varimp_cond

       pers_o        pers_c        pers_e        pers_a        pers_n         women     academics 
 1.656467e-02  2.611214e-03  2.404216e-03 -6.376514e-05  2.583184e-03  8.034409e-03  1.517586e-01 
          afd hospital_beds  tourism_beds           gdp      manufact       airport           age 
 1.486145e-01  1.175789e-02  2.029922e-03  9.026232e-02  1.075169e-02  1.350605e-02  2.393327e-01 
      popdens    onset_prev    slope_prev 
 4.531639e-02  3.076339e-03  1.778331e-02 

crf_meandiff_socdist_varimp_cond %>% as.data.frame() %>% 
  rownames_to_column('variable') %>%
  ggplot(aes(x=variable, y=.)) +
  geom_bar(stat = 'identity') +
  theme(axis.text.x = element_text(angle = 90))

Export data

ger_list_results <- list(cox_onset_prev, cox_onset_prev_ctrl, 
     lm_slope_prev, lm_slope_prev_ctrl, 
     cox_cpt_socdist, cox_cpt_socdist_ctrl, cox_cpt_socdist_ctrl2,
     lm_meandiff_socdist, lm_meandiff_socdist_ctrl, lm_meandiff_socdist_ctrl2)

results_names <- list('cox_onset_prev', 'cox_onset_prev_ctrl', 
     'lm_slope_prev', 'lm_slope_prev_coef', 
     'cox_cpt_socdist', 'cox_cpt_socdist_ctrl', 'cox_cpt_socdist_ctrl2',
     'lm_meandiff_socdist', 'lm_meandiff_socdist_ctrl', 'lm_meandiff_socdist_ctrl2')

names(ger_list_results) <- results_names

save(ger_list_results, file="ger_list_results.RData")

write_csv(df_ger_slope_prev, 'df_ger_slope_prev.csv')
write_csv(df_ger_cpt_socdist, 'df_ger_cpt_socdist.csv')